Plunge waveforms from inspiralling binary black holes 
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We study the coalescence of non-spinning binary black holes from near the innermost stable circular 
orbit down to the final single rotating black hole. We use a technique that combines the full 
numerical approach to solve Einstein equations, applied in the truly non-linear regime, and linearized 
perturbation theory around the final distorted single black hole at later times. We compute the 
plunge waveforms which present a non negligible signal lasting for t ~ lOOM showing early non- 
linear ringing, and we obtain estimates for the total gravitational energy and angular momentum 
radiated. 
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The next few years will mark the birth of a new field, 
gravitational wave astronomy. New extremely sensitive 
gravitational wave interferometers (LIGO and GEO) are 
nearing completion and should begin taking scientific 
data in about a year to be joined later by VIRGO and the 
LISA space mission. The expectation of very strong grav- 
itational wave emissions from the merger of black hole ~ 
black hole binary systems, and some indirect indications 
that these systems may be commonly generated in glob- 
ular clusters [Q makes them one of the most promising 
candidates for early observation. 

The interpretation of merger events, and in some cases 
even their detection, requires a theoretical understanding 
of the gravitational waveforms based on general relativ- 
ity. Several theoretical approaches have been developed 
for treating these systems. So far the post-Newtonian 
(PN) approximation, has provided a good understanding 
of the early slow adiabatic inspiral phase of these sys- 
tems, and extending PN calculations to the innermost 
stable circular orbit (ISCO) may be possible using re- 
summation techniques In its final moments, though, 
after the black holes are closer than ISCO, the orbital 
dynamics is expected to be replaced by a rapid plunge 
and coalescence. An understanding of the final burst of 
radiation coming from the black hole merger can signif- 
icantly enhance the detectability of systems with total 
mass AI > SSA/q |^ and even qualitative information 
is very important for developing search strategies It 
is generally expected that the dynamics near the ISCO 
can only be treated by a fully non-linear simulation of 
Einstein's equations. 

Thus far the numerical treatment of black hole sys- 
tems has proved difficult. A key limiting factor is the 
achievable evolution time after which the code fails due 
to numerical problems associated with the black hole sin- 
gularities. The first fully 3D simulation of spinning and 
moving black holes was performed in |^ for a 'grazing 
collision', where the black holes start out well within the 
ISCO. Recently, evolution times of about 9M-30M have 



been achieved However, these simulations, even be- 
ginning late in the plunge near the onset of quasi-normal 
ringing still generate only about two wave cycles. On 
the other hand, a third approach, the 'close limit' ap- 
proximation treating the late-time ring-down dynamics 
of the system as a linear perturbation of a single station- 
ary black hole has shown great effectiveness in estimating 
the radiative dynamics once the strong non-linear inter- 
actions have passed, for example in the case of grazing 
collisions without spin Q . 

In order to provide expectant observers with some es- 
timate of the full merger waveforms from binary black 
hole systems 'within a factor of two', and to guide fu- 
ture, more advanced numerical simulations, we have im- 
plemented an interface between full non-linear numeri- 
cal simulations and close limit approximations. This in- 
terface allows us to apply the numerical and close limit 
treatments in sequence, moving back the finite time in- 
terval of full non-linear numerical evolution to cover the 
earlier part of the plunge and then computing the com- 
plete black hole ring-down and the propagation of radia- 
tion into the wave zone with a close limit treatment based 
on the Teukolsky equation |^ . In we have developed 
the basic idea for the required interface and applied it to 
a model case, the head-on collisions of two black holes 
producing complete waveforms for the first time in full 
3D numerical relativity. 

In this letter we present the first theoretical predic- 
tions from non-axisymmetric binary black hole mergers 
starting from an estimate of the innermost stable circular 
orbit, based on a generalization of ||l0|] to include angu- 
lar momentum. We explicitly derive an astrophysically 
plausible estimate for the radiation waveforms and en- 
ergy which can be expected from a system of equal mass 
black holes with no intrinsic spin. We also estimate the 
duration of the plunge phase. 

Currently no genuinely astrophysical description of 
ISCO initial data for numerical simulations exists. As a 
reasonable starting point we will use approximate ISCO 
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data based on the effective potential method of as 
derived in for the puncture construction of black 
hole initial data. The solution to the general relativistic 
constraint equations is constructed within the Bowen- 
York, conformally flat, longitudinal ansatz. To locate 
the ISCO, the minimum in the binding energy along se- 
quences of constant (apparent) horizon area is studied 
as a function of the total angular momentum of the sys- 
tem. For puncture data the ISCO is characterized by the 
parameters 

TO = 0.45M, L = 4.9M, P = 0.335M, J = O.TTM^, (1) 

where to is the mass of each of the single black holes, 
M is the total ADM mass of the binary system, L is 
the proper distance between the apparent horizons, P is 
the magnitude of the linear momenta (equal but opposite 
and perpendicular to the line connecting the holes), and 
J is the total angular momentum. 

Our full numerical evolutions are carried out using the 
standard ADM decomposition of the Einstein equations. 
The lapse is determined by the maximal slicing condi- 
tion, and the shift is set to zero. Our biggest runs had 
512^ X 256 grid points and a central resolution of M/24. 
We use the Cactus Computational Toolkit to implement 
these simulations . Our typical evolution times reach 
T K, 15Af. After this time, the 'grid stretching' associ- 
ated with the singularity avoiding property of maximal 
slicing crashes the code. This is consistent with the evo- 
lution time for grazing collisions of about 30M, because 
for TSCO' runs the time scale for grid stretching is de- 
termined by the individual black holes of mass 0.45M 
that are present for most of the run, and not by the fi- 
nal black hole of mass \M that quickly appears in the 
grazing collisions. 

The fundamental problem in determining initial data 
for the close limit perturbative approach is to define a 
background Kerr metric on the later part of the numer- 
ically computed spacetime. We look for a hypersurface 
and a spatial coordinate transformation such that the 
transformed numerical data on this hypersurface is close 
to a Kerr metric in Boyer-Lindquist coordinates to which 
all perturbative quantities refer to. The first order gauge 
and tetrad invariance of the perturbative formalism im- 
plies that the results will not depend strongly on small 
variations in our procedure. 

We start with the estimate that the background Kerr 
black hole is given by the parameters M and a of the ini- 
tial data. Then we correct those values with the informa- 
tion about the energy and angular momentum radiated. 
This procedure quickly converges. The slice is defined 
by identifying the numerical time coordinate with the 
Boyer-Lindquist time, since we have found that in the 
late stages of the numerical evolution the maximal slic- 
ing lapse approaches quite closely the Boyer-Lindquist 
lapse {—g^Kerr) outside the black hole. To obtain the 



r coordinate we compute the invariant X = Cap-ysC"^^ , 
where Cap^s = CafS^s + {i/2)eaf3p\C''^g, is the self-dual 
part of the Weyl tensor. Since in Boyer-Lindquist coor- 
dinates X = 3M^/ (r — ia cos 0)^ for the Kerr background 
we can invert this relationship to assign a value of r to 
each point of the numerically generated spacetime. For 
the coordinate 9 it turned out to be sufficient to adopt 
the corresponding numerical coordinate. Since the bi- 
nary system carries a non-negligible amount of angular 
momentum J producing frame dragging effects, we sup- 
plement the Cartesian definition of (j) with a correction 
that makes the gr,p component of the metric to vanish 
(primes denote numerical coordinates) 

(l> = (t>' + j{9'r'^'l9'^'^')dr'. (2) 

The Teukolsky equation is a complex linear wave equa- 
tion for the two degrees of freedom of the gravitational 
field in a Kerr background. It requires as initial data 
the Weyl scalar -04 = —Cap-ysn^fh^rfifi^, with the usual 
notation for the null tetrad, and its time derivative dtipi 
(see for the ADM decomposition of these quantities). 
First, a tetrad is constructed from the basis vectors of the 
numerical coordinates, and a Gram-Schmidt procedure is 
used to ensure orthonormalization with respect to the full 
numerical metric. The tetrad is then rotated to approx- 
imate the required background tetrad. This procedure 
determines ip4 as a linear combination of all the other 
five Weyl scalars ■04, 3, 2,1,0 with coefficients depending on 
the background coordinates r and 9, and on M and a 
JlGf . With the Cauchy data 04 and dttp4 extracted in 
appropriate coordinates, we can proceed with the evolu- 
tion by numerically integrating the Teukolsky equation 
as described in Here one can implement all the de- 
sired features for stable evolution, excision of the event 
horizon, mesh refinement through the use of the tortoise 
coordinate r*, non- vanishing background shift, imposi- 
tion of consistent boundary conditions on 04, etc. The 
perturbative description is then able to efficiently follow 
the evolution of the system forever. 

A key question is at what time we can actually make 
the transition from full numerical to perturbative evo- 
lution. One of the tests we perform is to compute 
at every IM of evolution the speciality index S = 
27J^^/T^ [Q, a combination of curvature invariants 
(J = C'apysC^lx^''^"'^)- For the algebraically special 
Kerr spacetime S — 1, and the size of deviation from 1 
provides a guide on how close a numerical spacetime is to 
Kerr. A complementary experiment is to do a full numer- 
ical evolution, during which we extract the data for the 
Teukolsky equation every IM. Computing the Teukol- 
sky evolution starting at different times, we can check 
whether the final results (radiated energies and wave- 
forms) depend on the time T at which the transition took 
place. If the binary system has reached a regime where all 



2 



ISCO, r*=31M, z-axis, m=+2 



P-Sequence 



15 



10 



: Common Apparent Horizon: : 



0.2 0.4 0.6 0.! 

P/P 




FIG. 1. Linearization time and total radiated energy. 

further evolution can be described by the linearized Ein- 
stein equations, the final results should be independent 
of the choice of the transition time T. This constitutes 
an important built-in self consistency test of our method. 

As mentioned above, we have successfully tested the 
basic procedure for the head-on collision of black holes 
pO| . In we report on non-trivial consistency checks, 
e.g. quadratic convergence to vanishing gravitational ra- 
diation, which our method passes for the evolution of 
Kerr initial data. Now, in order to better understand the 
physics of the plunge we have designed a set of sequences 
approaching the ISCO by varying one of its physical pa- 
rameters. Several different sequences are possible, which 
we will discuss in a longer paper. Here we report only on 
the 'P-sequence' for which we keep the separation con- 
stant aX L = Ljsco = 4.9M, but vary the linear momen- 
tum, P/ Pisco ~ 0, |, |, |, 1, y|. With this sequence we 
can vary continuously from the head-on collisions treated 
effectively in to the case we are most interested in, 
P — Pisco- None of the elements of this sequence starts 
in the close limit regime. Referring to Fig. |^, we ob- 
serve that the minimal time of full numerical evolution T 
needed to switch to perturbation theory initially grows 
roughly linearly with increasing momentum, with a fur- 
ther deviation towards longer times when approaching 
the ISCO. We also show the time at which a common ap- 
parent horizon appears, if it actually appears during the 
achievable 15M of full numerical evolution. In our sim- 
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FIG. 2. ISCO waveform for two transition times compared 
against the result with no full numerical evolution, T — 0. 



ulations this provides an upper limit to the linearization 
time. We note that with the current method P = Pisco 
is a marginal case. Based on the S invariant and the tran- 
sition time experiment, linearization for the ISCO hap- 
pens in the range T = 11M-15M. For P/ Pisco = i| 
the S invariant does not indicate linearization much be- 
fore T = 15M. As shown in Fig. ||, the total radiated 
gravitational energy grows quadratically with the mo- 
mentum P for small values, as one would expect from 
dimensional arguments and from extrapolation of close 
limit results The error bars are based on variations 
for different transition times. 

Fig. ^ shows the real part of 1/^4 as seen by an observer 
located at radial coordinate r* = 31M along the orbital 
pole. We display the leading m = 2 mode, since we 
decompose the Weyl scalar into e*™*^ modes. The wave- 
forms are obtained for 0, 10, and IIM of full numerical 
evolution plus 120A/ of linear evolution. The T = Q 
waveform is displayed only to illustrate the importance 
of the full numerical evolution. As is suggested by com- 
paring the waveforms for T = lOM and IIM, we have 
precisely determined the first part of the waveform last- 
ing up to about t = BOM. Beyond that the agreement 
is less precise. In the final region some details of our re- 
sult, in particular the phase, are still sensitive to elements 
of the approximation procedure. Nevertheless, we judge 
that we have met our goal of a waveform estimate within 
a 'factor of two' even in that region. The exact shape 
of the waveforms after t — 60M may still change when 
more sophisticated numerical techniques are applied. 

Let us summarize some of the robust features of our 
waveform. The duration of the most significant part of 
the waveform is roughly lOOAf , lasting for about twice 
as many wave cycles as in the head-on collision case. 
There is a 'ring-up' with increasing amplitude until about 
t « 60M followed by a 'ring-down'. If we define the dura- 
tion of the plunge by the time interval from the beginning 
of the signal to the onset of the ring-down in Fig. |^, we 
obtain a plunge time of 30 Af. This is comparable to 
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the plunge time derived by different means in PN calcu- 
lations in the range from AIM to 77 M ||]. For the type 
of initial data we use the plunge time is estimated to be 
equal to the orbital period at the ISCO in [ pl]jT^ ], which 
is 37M and again consistent with what we find. 

For an orbital period of 37M we expect to see 
quadrupole radiation at half that period. In fact, a fre- 
quency decomposition of the waveform shows a dominat- 
ing component with periods close to 12M and 19M, the 
two most weakly damped quasi-normal oscillation fre- 
quencies for m = ±2 [|l9| of our final black hole with 
a = 0.8M. This suggests that at around 19M there is a 
mixture of quasi-normal ringing and orbital components. 
For a system with total mass 35Mq our two principal fre- 
quency components correspond to frequencies of roughly 
300Hz and 4:75Hz, which are within the sensitive range 
of typical interferometric gravitational wave detectors. 
Note that there is a discrepancy of roughly a factor of 
two in the orbital period at the ISCO between the initial 
data analysis and PN calculations, i.e. our initial data 
corresponds to tighter configurations than are found in 
PN studies. 

Let us point out that the observed time for lineariza- 
tion of less than 15 A/ is significantly shorter than the 
« 30Af it takes until ring-down. The onset of linear evo- 
lution need not immediately result in a ring-down wave- 
form. The linearization time is also less than the time un- 
til a common apparent horizon appears, but recall what 
sets the range of validity of the close limit approxima- 
tion. It is not the presence of a common horizon, but 
rather that the black holes sit in a common gravitational 
well, which occurs earlier. The briefness of the non-linear 
phase of the plunge is good news, because this is the tech- 
nical reason why we can perform these simulations with 
the current techniques. Still, we want to emphasize that 
numerical relativity was essential for achieving these re- 
sults. For T = 11 M, the first wave cycle is determined 
precisely to within 1% of a wave cycle, while the T — OM 
waveform with no full numerical evolution has a very dif- 
ferent appearance and is roughly 90 degrees out of phase. 

We estimate the total radiated energy after ISCO to be 
w 3% of the total ADM mass coming almost entirely from 
the m = ±2 modes. This is larger than the 1.4% obtained 
by extrapolating PN results and the 1-2% obtained by 
extrapolation of the close limit jsj . The radiated angular 
momentum is a delicate quantity to compute involving 
correlations of waveforms [po| , we estimate the angular 
momentum loss to be around 0.2%. This confirms the 
expectation that not much angular momentum is lost 
during the plunge and ring-down. 

In conclusion, our approach makes it for the the first 
time possible to study the fundamentally non-linear pro- 
cesses taking place during the final plunge phase of the 
collision of two well-separated black holes, starting from 
an estimate for the ISCO location. The results presented 
here show that full 3D numerical evolution is essential to 



describe the non-linear interaction of binary black holes. 
The interface to a linear evolution of Einstein equations 
allows us to target the full numerical evolution where it 
really matters. 

Ours are the first but certainly not the final numeri- 
cal results for black hole systems near the ISCO. Further 
work will be done to extend the duration of the numerical 
simulations and to move toward more astrophysically re- 
alistic initial data descriptions, e.g. given by an interface 
to the post-Newtonian approximation. The technology 
we have developed makes it possible to explicitly study 
the physical appropriateness of various initial data de- 
scriptions, both by comparison studies, and by moving 
to more separated 'pre-ISCO' configurations. As another 
astrophysical application of numerical relativity one can 
now begin to characterize the effects of black hole spins 
and relative mass differences on the radiation waveforms. 
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